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We have used numerical modeling and a semi-analytical calculation 
method to find the low frequency value S/(0) of the spectral density of fluc- 
tuations of current through ID arrays of small tunnel junctions, using the 
"orthodox theory" of single-electron tunneling. In all three array types stud- 
ied, at low temperature (fcgT <C eV), increasing current induces a crossover 
from the Schottky value 5/(0) = 2e(I) to the "reduced Schottky value" 
Si(0) = 2e{I)/N (where ./V is the array length) at some crossover current 
I c . In uniform arrays over a ground plane, I c is proportional to exp(-XN), 
where A -1 is the single-electron soliton length. In arrays without a ground 
plane, I c decreases slowly with both N and A . Finally, we have calculated 
the statistics of I c for ensembles of arrays with random background charges. 
The standard deviation of I c from the ensemble average (I c ) is quite large, 
typically between 0.5 and 0.7 of (I c ), while the dependence of {I c ) on N or A 
is so weak that it is hidden within the random fluctuations of the crossover 
current. 
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I. INTRODUCTION 



Single-electron tunneling (for general reviews, see Refs. |l]fj) is one of the most active 
areas of solid state physics research, but it has several key problems that have not yet been 
addressed in detail. One of these problems is the dilemma of the discreteness/continuity of 
electric charge transfer^. Several single-electron tunneling phenomena can be understood as 
an interplay between discrete transfer of electric charge via electron tunneling, and continu- 
ous transfer in ordinary diffusive conductors. For example, a tunnel junction biased by a dc 
current may generate single-electron tunneling (SET) oscillations with average frequency! 

This effect may be interpreted as a gradual accumulation of continuous charge on the junction 
capacitance, followed by the abrupt passage of one electron through the junction, as soon 
as the accumulated charge has reached a threshold level! Q t = ±e / 2 . 

However, if the charge transfer in the external circuit (fixing the current /) is discrete, 
SET oscillations do not exist. This can be seen from the following general formula for the 
SET oscillation linewidth0H: 

T s =(^yS I (0) (T s <f s ), (2) 

where Sj(f) is the spectral density of the bias current fluctuations. For example, if the 
current is fixed using another tunnel junction, then at low temperatures the fluctuations 
obey the Schottky formula 

Sj(0) = 2e(I) , (3) 

and Ts > fs, that is, SET oscillations are completely smeared by the current fluctuationsi. 
On the other hand, in macroscopic diffusive conductors, the current noise may be much 
lower than the Schottky value, and SET oscillations (and very similar "Bloch" oscillation 
may have a relatively narrow bandwidth - see, e.g., experimental measurements in Ref. [5[ 
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More generally, discreteness of charge transfer will certainly be one of the central is- 
sues facing the emerging nanoelectronics. This is why it is important to formulate the 
conditions under which the transport of charge through a conductor may be considered as 
(quasi)continuous, i.e., having discreteness 5Q <C e. The definition of 5Q can be most readily 
introduced in the most interesting case of negligible thermal and quantum fluctuations: 

k B T < [e(V), eV t ] , (4) 
G<£, (5) 

where V t is the Coulomb blockade threshold of the conductor and G its effective conductance. 
However, even in this simplest case, the definition depends essentially on the characteristic 
time scale r of an experiment: 

A. If fsT 3> 1 (i.e. either the time scale r is large, or the dc current is high, or both), 
SQ may be defined as follows: 

5_Q Sj{0)_ 

e 2e(I) ' 1 ' 

In fact, if the charge Q transferred through a system may be presented as a Poissonian 
series of jumps of fixed height 5Q, then repeating the well-known derivation of the Schottky 
formula we arrive at Eqn. (||). If the jump height is random as well, Eqn. (|6|) is still 
applicable as an estimate of the average jump height. 

B. In the opposite limit, when r is much shorter than the average spacing between the 
charge jumps (fg 1 ), we are essentially dealing with the Coulomb blockade regime. In this 
case an adequate definition of 5Q is as follows: 

SQ = C m V t , (7) 

where C{ n is the effective input capacitance of the system of interest. (5Q given by this 
formula is the fraction of the initial electric charge of the system which cannot relax, due to 
the Coulomb blockade; for a diffusive conductor 5Q — > 0, while for a single tunnel junction 
5Q = e/2.) 
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One of the most interesting systems capable of quasicontinuous charge transfer is the ID 
array of small tunnel junctions — see, e.g., the review in Ref. |6|. The key property of such 
an arrayi is that each additional electron inserted into one of its islands creates a series of 
gradually decreasing polarization charges, and hence may be considered as a "single-electron 
soliton" with a characteristic size M which may be much larger than one island. As a result, 
when an electron is drifting along the array, tunneling between the neighboring islands, the 
effective charge Q transferred through the external electrodes — which is essentially what 
interests us — changes in jumps of scale 5Q ~ e/M which may be much smaller than e. 

This does not mean, however, that the condition M > 1 is sufficient for the quasicon- 
tinuous transfer of charge in the arrays. In fact, the calculation of the Coulomb blockade 
threshold for arrays with uniform! and randomHill background charges have shown that, 
under definition B, charge transport in uniform arrays with vanishing background charge 
cannotbe considered as continuous (5Q ~ e). The reason is that a relatively strong Coulomb 
blockage results from the pinning of single-electron solitons by the sharp edges of the ar- 
ray. Presently, we know only two cases when charge transport through a uniform array is 
quasicontinuous (5Q <C e) according to the definition B: 

1) An array with arbitrary capacitances, but with the background charge of all the islands 
equal to ±e/2. (In this case the Coulomb blockade threshold vanishes, and 6Q — ► 0.) 

2) An array with 1 <^ N <^ 2M and random background charges. (In this case00 
V t ~ O.heN 1 / 2 / C , where C is the tunnel junction capacitance, while Cj n = C/N, and hence 
5Q/e « 1/2N 1 ' 2 < 1.) 

However, even in these cases, the charge transfer is not automatically continuous in the 
sense of criterion A. If the Coulomb blockade is finite, and the array is driven with dc voltage 
V slightly above the blockade threshold V t , one of the junctions presents a bottleneck to the 
single-electron soliton drift along the array. As a result, the passage of an electron consists 
of a long wait at the bottleneck junction, followed by a rapid burst of transitions through 
the remaining junctions of the array. At V — » Vt the statistics of these bursts is always 
Poissonian, and the charge transferred by each burst is equal to e, so that the shot noise is 



well described by the Schottky formula. 

Thus, current noise in ID arrays presents an important problem. To our knowledge, this 
problem has previously only been solved^ for a very particular case of uniform arrays with 
M — > oo (zero stray capacitance) and zero background charge. The objective of this work 
was to calculate Sj(0) (and hence the effective discreteness of charge transfer for long time 
intervals) for a much broader range of array parameters. 



II. BASIC FORMULAS 

We have considered arrays consisting of iV — 1 small metallic islands connected by N 
tunnel junctions, and flanked at either end by dc voltage-biased electrodes (Fig. [I]). Under 
the conditions expressed by Eqs. (f|), (J5|), we can ignore the effects of cotunneling and 
of thermally activated tunneling. Current flow in single-electron arrays may be analyzed 
in terms of the junction tunneling rates, T^. According to the orthodox theoryi, at zero 
temperature, 

«A^,>0, 

r 4 ,= e2 (8) 

I AWij < , 

where AWy is the drop in the free (electrostatic) energy, caused by the tunneling event. 
The drop in energy due to an electron tunneling from island % to island j can be written 



e 2 CZ x + e 2 C~ x 

Y 



AWv = e 2 C^ ZL + e[(j)j _ ^ , (9) 



where <fi is the vector of the electrostatic potential of the islands before the jump, and the 
matrix C _1 of inverse capacitances is defined by the following equation: 

& = E c vfa + & > * = E ■ ( 10 ) 

j&sl kdext 

Here the matrix C represents capacitances between islands and external terminals with 
potentials Vk- If these potentials do not change in time, the probability that the system 
preserves its charge state can be expressed explicitly: 



P(*) = exp(-r(*-* )), (11) 

where t Q is the time of the preceding tunneling event, while T is the total rate for all possible 
tunneling events: 

r= £ *V ( 12 ) 

{ij}ejrf 

III. NOISE COMPUTATION 

The preceding relations were incorporated into our main computational tool, a C++ 
program called msol that uses a Monte Carlo algorithm! to simulate the flow of current in 
dc voltage biased ID arrays. The basic unit of calculation in mso is the "current run", in 
which charge flows through the array until a user-specified total charge Q is transferred. To 
calculate S/(0), mso loops through a user-specified number N r of current runs, each starting 
with the same initial charge state and ending when the total charge transferred equals Q. 
The same random number generators^! for time and jump location are used continuously 
through all loops. To the extent that the random numbers constitute a "good" quasi-random 
series, the ensemble of current runs represents a statistical ensemble of independent, initially 
identical systems. 

We may calculate S/(0) from the statistical properties of the time Tq taken by each run, 



asL 2 



Sr(0) ^Q Var(r q ) 

2e(I) ~ e (Tq) 2 ' 1 ) 

where I = Q/Tq, and the angle brackets and variance refer to our statistical ensemble of 
current runs. Equation ([13]) is exact only for infinite Q and N r ; since the jumps are not 
completely independent, the accuracy of this formula should be determined experimentally. 
Figure |] shows a typical dependence of S/(0)/2e(7) on Q and N r . The results for Q > lOOOe, 
iV r > 1000 seem to be accurate to within 10% of the asymptotic value. In this paper, we 
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used the parameters Q = lOOOe, N r = 1000 for calculating shot noise in arrays without 
background charge, and assigned 10% error bars to these numbers. 

For arrays with random background charge, each point was calculated for 50 different 
realizations of the background charge for each circuit, using the parameters Q = 200e and 
N r = 200, to keep the simulation time within reasonable limits. Although these calculations 
are therefore less accurate, perhaps only to within 20% of their asymptotic value, this 
inaccuracy was overshadowed by the overall spread in shot noise values among the different 
background charge realizations. 

Calculation (CPU) times in mso scale as ~ N a QN r , with a slightly larger than 2. A 
typical calculation with N = 20, Q = lOOOe, and N r = 1000 takes around 400 seconds 
of CPU time on an AlphaStation 250 (266 MHz Alpha, Digital Unix 4.0b, Digital cxx) or 
around 950 seconds on a Linux PC (120 MHz Pentium, RedHat Linux 2.0.30, Gnu C++). 



IV. CROSSOVER CURRENT 

Varying the bias voltage across an array, we have calculated the average current and spec- 
tral density as functions of applied voltage and have made parametric plots of 5/(0)/2e(7) 
vs. (7) (Fig. |H). We will refer to these plots as 5 — I curves. 

The most immediate, universal result of our calculations is the crossover of 5/(0)/2e(7) 
from 1 to 1/N with increasing current. This result can be understood as follows. As argued 
in section |, 5/(0) near threshold is dominated by the Poissonian statistics of tunneling 
through a single bottleneck junction, and is thus given by the Schottky formula: 

5,(0) 



-1. (14) 



2e(I) 

At high voltages, however, a large number of charge states becomes available for tunneling 
through each junction. Though the tunneling rate for each of these states may be affected 
by the state of neighboring junctions, these effects are averaged out, since the voltage depen- 
dence of the rate of tunneling through each junction is linear at AWij > (see Eqn. (^)). 
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Under these conditions, current noise through each junction is described by the Schottky 
formula, S?\o) = 2e(7) . S ince we may transform current noise into voltage noise by the 
square of the dynamic resistance (which at V — > oo just equals R), we write the voltage 
noise of a single junction as 

S { y\0) = 2e(I)R 2 . 

The total voltage noise Sy is the simple sum of the noise of the individual junctions Sy\ 
while the total current noise is finally obtained from Sy via the total array resistance NR: 

Sy(0) 1^ = 2e{I)R 2 N, (15) 



(For the particular case N = 2 this equation has been discussed in Ref. TA.) Thus the 
crossover in Sj(0)/2e{I) from 1 to 1/N with increasing current could be expected; what 
was really surprising for us is that in most cases this crossover takes place very close to the 
Coulomb blockade threshold, where the array I — V curve is still not quite linear, and hence 
the arguments given above cannot be taken too seriously. 

In order to describe the crossover quantitatively, we may define the crossover current I c 
as a value at which Sj(0)/2e(I) is midway between these two limits, on a logarithmic scale: 

(17) 

We have written a Perl script called si curve to automate the extraction of the crossover 
current I c from the S — I curves. While invoking mso, sicurve continuously adjusts the 
change in bias voltage between successive points, in an attempt to produce a series of evenly 
spaced points in the S/(0)/2e(7) vs. (7) plane (see Fig. §). This is an important practical 
technique for generating S — I curves on circuits with random background charge, since the 
relationship between voltage, current, and spectral density can be quite irregular. 
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V. MODEL 1: ARRAYS NEAR GROUND PLANE 

Our first case was the simplest mo del of a uniform, symmetrically biased array near a 
ground plane with no background charges (Fig. lb). The direct capacitance matrix in this 
model is tridiagonal, and is described by one dimensionless parameter, the ratio Cq/C, where 
C is the junction capacitance, and Cq is the "stray" capacitance between an island and the 
ground plane. In this model, the reciprocal length scale A = 1/M of the single-electron 
soliton is determined as 

A = cosh" 1 (C /2C-l), (18) 



and in the most interesting limit of Co <C C, M = y^C/Co 3> 1. 

Surprisingly, our numerical results (Fig. f|) show that for all values of A and N, data for 
I c fall roughly on a single universal curve, with I c oc exp(— NX/3) for A > 10. For a fixed 
product NX, there is a relatively weak decrease of I c with increasing N. 

To understand this unexpected result, we began looking for an analytic expression for 
current noise for the case when the passage of charge through the array consists of a fixed, 
repeated sequence of tunneling events. If such a sequence is repeated n = Q/e times, we can 
write the average total time as a sum of average times for each jump in the fixed sequence: 

- n - 

(T Q ) = n^iSU) = J^i 1 ' ( 19 ) 

i=l i=\ 

where T is the total rate at each stage of the process. Since the time of each jump follows 
Poissonian statistics, 

Var(^) = Yf , (20) 
and since the jump times are independent of one another, 

n - o N 

Var(T Q ) = ^J2 Var (^) = 7 E T * 2 > ( 21 ) 

i=l i=l 



we arrive 



ive at a simple formula0: 
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TV 

-2 



Si(0) 
2e(J) " 



E^ 

(22) 



E^ 1 

v 1=1 / 

This formula was obtained earlierB for a particular case Co = 0. In particular, it shows 
that the spectral density can be dominated by bottleneck points, where the rate Tj is much 
lower than average. 

In order to find these bottleneck points for the present case of arrays near a ground 
plane, let us examine the energy profile created by these arrays for tunneling charges. The 
potential created by an electron in an array with externals at both ends isEl 

e / _| W |A sinh[(iV - n')\] + e^ N ~^ sinh(nA) | 

C^l sinh(iVA) J ' {l6) 

where n e {0,1... iV} is the position of the electron in the array, n' is the measurement 

position, and C e // = a/ Co + 4CC . For our symmetric bias (V 2 = — Vi = V/2), the potential 

created by the external electrodes iM 

V^ H{ N/2-n)X] 
r V ; 2 sinh(ATA/2) V ; 

Numerical simulations show that in symmetrically biased arrays with A > 1, the basic 
tunneling scenario near threshold is the passage of electron-hole pairs. The components of 
the pair enter at opposite ends of the array, move towards each other, then annihilate near 
the center (see Table 1). We can write the energy of the electron-hole pair as 

W(ni, n 2 ) = W (ni) - e<j) e (nx) + W {n 2 ) + e0 e (n 2 ) - e<p s (nx, n 2 ) , (25) 



where rii and n 2 are the positions of the electron and hole, respectively, and 

cosh[(JV-2n)A] -e _iVA ' 



Wo(n) = 



sinh(iVA) 



(26) 



2 2C eff 
is the self-energy of an electron or hole. 

Figure |] shows the results of calculation of the crossover current using Eqs. (|25|) and 
(f2~2|) in a fixed scenario picture, where the electron and hole enter the array one right after 
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the other, and then take turns tunneling towards each other in a symmetric manner. For 
iVA > 15, the results of this semi- analytical calculation match the Monte Carlo results very 
closely. (For NX < 15, a difference appears, increasing with N.) This means that we can 
analyze our problem, at least approximately, by examining the energy profile in our fixed 
scenario. A straightforward analysis of Eq. (ESj) shows that, for NX ^> 1, there is a slow 
point (a minimum in AWij) when both electron and hole are ~A r /3 islands from their 
respective edge of the array. Here both members of the electron-hole pair are far from the 
edges and cannot be pushed strongly by the external voltage, yet are not close enough to 
attract each other strongly, either. At this point, T scales as exp(— NX/3). According to 
Eq. (p2|) , this leads to a similar dependence of J c , at least in the limit iVA ^> 1. 

Some difference between the fixed-scenario calculations and the Monte Carlo simulations 
can be readily explained by the observation that frequently the tunneling process is some- 
what more complicated than the exactly alternating electron and hole motion sequence — see 
Table 1. Figure |5| shows a typical pattern of the total tunneling rates. 

These figures show that sometimes there are three (rather then one) slow points with 
low T. The first slow point is before any charge enters the array (see labels (a) in Fig. |5|). 
After the first charge enters the array, the external voltage pushes it quickly away from the 
edge. There is then typically a second slow point when that charge reaches a distance of 
roughly iV/3 islands from the external it emerged from (see labels (b) in Fig. |5|). Then the 
opposite charge enters the array from the other side, also hopping quickly away from the 
edge. The third slow point usually follows when the opposite charge is roughly N/3 islands 
away from its external (see labels (c) in Fig. |5|). Only the first two points are taken into 
account in our fixed scenario. Since the second slow point is frequently important, the good 
agreement between our Monte Carlo and fixed scenario calculations is somewhat puzzling. 
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VI. MODEL 2: ARRAYS WITHOUT GROUND PLANE 



Here we examine the case of an array with islands between semi-infinite external elec- 
trodes (Fig. ||a), without a ground plane. In the absence of a ground plane, mutual capaci- 
tances other than the junction capacitances may become important. 

To model the electrostatics of such arrays, we have used two methods. In the first 



use it to calculate the full capacitance matrix for the array numerically using FastCapB. 
Since the electrostatics is rather insensitive to the exact shape of the islandsS, they may be 
modeled by cubes of side length a. The resulting capacitance matrix (an example is given 
in Table |T|) was used for the Monte Carlo simulation of noise, as described above. 

In the second method, we used a simple heuristic approximation for the single-electron 
soliton potential at distance m = \n — n'\ in a long array, found in Ref. [17], 



with a ~ k ~ 1. This formula describes a crossover from an exponential decay at short 
distance (Am < 1) to a Coulomb-law 1/r decay at large distance. The effect of the external 
electrodes was described by the usual image charge method; in our case, with two electrodes, 
it involves an infinite series of images. As a result, the full single electron potential, 0(?(n, n'), 
may be expressed as the sum of an infinite series over all image charge contributions, and 
the self energy can be written as 





(27) 



Wo h 



(n) 



e<f)^(n, n) 
2 



[W a (2N(i - 1) + 2n) + W a (2Ni - 2n)} 
2 



On the other hand, the form of the external potential is quite simple, 



(j) h e {n) ~ eV(N -2n)/2N . 



The energy for an electron-hole pair can then be written just as in Eqn. This energy 
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was used for the fixed-scenario calculation of noise, similar to that described in the previous 
section. 

Results from both methods are shown in Fig. |7|. The Monte Carlo simulation shows 
that, unlike in arrays near a ground plane, I c does not follow a universal dependence on 
iVA. Rather, the I c values follow a common curve depending mostly on A alone, decreasing 
weakly with the array length N. Even the A dependence is weak, compared to the results 
for arrays near ground plane: for iV = 20 in the range of A from 0.25 to 1, I c drops by a 
factor of ~ 3 in arrays with no ground plane, whereas it drops by more than two decades in 
arrays with ground plane (Fig. ^). 

The fixed scenario results for J c match the Monte Carlo results fairly well in terms of the 
shape of the A dependence of the curves. For iV = 10, the magnitude of the results are also 
in fairly close agreement. However, the fixed scenario results exhibit a stronger decline with 
N than the Monte Carlo results. For N = 60, the fixed scenario values of I c are between 3 
and 8 times smaller than the corresponding Monte Carlo values. 

Generally, it is easy to understand why the fixed scenario values for I c fall below the 
values calculated with Monte Carlo simulation: randomness of jump location, which is 
ignored in the fixed scenario calculation, can be thought of as an additional source of noise. 
As the general level of noise in the crossover region increases, so does the crossover current. 
However, we are still in need of a simple interpretation of the N and A dependence of the 
shot noise in arrays without ground plane. 

VII. MODEL 3: ARRAYS WITH RANDOM BACKGROUND CHARGE 

Returning to arrays near ground plane, with their simple electrostatics, we have explored 
the behavior of shot noise in the presence of random background charges on the islands, go- 
These charges can represent, for example, the effect of charged impurities in a substrate. 
The charge ^ placed on each island was randomly selectedi, using a uniform probability 
distribution between — e/2 and e/2 . (Any integer part of the background charge would 
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immediately be compensated by trapping one or a few tunneling electrons or holes.) 

The S — I curves which stem from the Monte Carlo simulation of such arrays still show 
the progression from shot noise to suppressed shot noise, but at larger values of A tend to 
feature strong, irregular peaks (Fig. |3]) where curves for arrays without background charge 
were smooth. Occasionally these peaks in the S — I curve cause it to cross the 1/yN line 
more than once. In such cases, we arbitrarily define I c as the lowest crossing value. It 
turns out that the variation in the general location of I c among different instantiations of go 
is large enough to make the distinctions among different crossing points on a single curve 
irrelevant. 

Figure |8] shows the average value of I c and its standard deviation of over ensembles of 
50 different % realizations. Strikingly, for A > 1, I c for all values of N appears to fall on 
the same, almost flat curve. This behavior is in sharp contrast to results for similar arrays 
without background charge (Fig. This quasi-universal value of I c is ~ 5 x 10~ 3 e/ RC. 
Although there appears to be weak downward trend in (I c ) with increasing A for A > 1.5, 
this trend is virtually hidden within the relatively large standard deviation. 

Let us try to comprehend this result. Simulation shows that in such arrays the typical 
tunneling process near the crossover is due to several (rather than one) electron-hole pairs 
moving simultaneously in the array. (Due to this multiplicity, a fixed-scenario calculation of 
noise would not make sense.) This effect is easy to explain. The electric potential induced 
by the background charges creates a series of charge traps. So, as the applied voltage is 
increased, the first several charges to enter the array are trapped, forming a "sandpile"! of 
charge (see Fig. |^ ). The final energy profile in the sandpile is still random, though the 
maximum possible energy change AWij is now upwardly bound by AW max ~ e2 /2C e //. 
The remaining disorder is, however, strong enough to overcome the interaction of distant 
charges, which is exponentially weak in arrays over a ground plane. 

This argument also explains why the crossover current is independent of the array length. 
We are, however, still in need of an analytical theory which would explain the virtual in- 
dependence of the average crossover on A , and also would predict its universal value cited 
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above. 

VIII. CONCLUSION 

We have used both Monte Carlo and fixed scenario techniques to calculate the low- 
frequency current noise for three different models of ID single-electron-tunneling arrays. 
Within each of the three models, we find a crossover of the spectral density of current 
fluctuations, S/(0), from the Schottky value 2e(J), to the "reduced Schottky" value 2e(I)/N, 
with increasing current. The crossover can be well characterized in terms of the crossover 
current J c , which may be said to mark the onset of quasicontinuous charge transfer. The 
particular behavior of I c as a function of N and A depends on the interaction of electrons 
within the array and on the interaction of electrons with the external field. 

For arrays near ground plane, with no background charge, the crossover current exhibits 
a universal behavior that is a function of only the product iVA , i.e. of the ratio of the array 
length N to the length M = A -1 of single-electron solitons. At iVA > 10, the dependence 
is exponential: I c oc exp(— NX/3). Our analysis has shown that this behavior is the result 
of the exponential decrease of the soliton interaction with the external electrodes and its 
counterpart in the electron-hole pair. 

In arrays without a ground plane, I c is almost independent of N and exhibits a nearly- 
universal weak decrease with A. We believe that the substantial difference in results between 
this model and the previous model is due to the long-range electrostatic interactions, which 
were screened by the ground plane in the previous model. 

Finally, in arrays with random q Q , we have found that the crossover takes place at a 
nearly universal value of current, I c ~ 5 x 10~ 3 e/ RC. The absence of noise dependence on 
the array length can be readily explained as a result of the random potential created by 
the background charges, which overwhelms long range order in the arrays. However, the 
independence of noise on the single-electron soliton length still has to be explained. 



15 



IX. ACKNOWLEDGMENTS 



The authors thank D.V. Averin and A.N. Korotkov for fruitful discussions, and A. Huq 
for technical assistance. 



1 D. V. Averin and K. K. Likharev, Mesoscopic Phenomena in Solids (Elsevier, Amsterdam, 1991), 
p. 173. 

2 Single Charge Tunneling, edited by H. Grabert and M. Devoret (Plenum, New York, 1992). 

3 D. V. Averin and K. K. Likharev, J. Low Temp. Phys. 62, 345 (1986). 
4 K. K. Likharev and A. B. Zorin, J. Low Temp. Phys. 59, 347 (1985). 
5 L. S. Kuzmin and D. Haviland, Phys. Rev. Lett. 67, 1161 (1991). 

6 P. Delsing, in Ref. 2, p. 249. 

7 K. A. Matsuoka, 1997. Software available at http : / /hana. physics . sunysb.edu / mso / . 



8 N. S. Bakhvalov, G. S. Kazacha, K. K. Likharev, and S. I. Serdyukova, Sov. Phys. JETP 68, 
581 (1989). 

9 A. A. Middleton and N. S. Wingreen, Phys. Rev. Lett. 71, 3198 (1993). 

10 J. A. Melsen, U. Hanke, H. O. Midler, and K. A. Chao, Phys. Rev. B 55, 10638 (1997). 

11 A. N. Korotkov, unpublished (1996). 

12 A. N. Korotkov, Phys. Rev. B 50, 17674 (1994). 

13 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C 
(Cambridge Univ. Press, Cambridge, 1992), pp. 278-282. 

14 M. J. M. de Jong and C. W. J. Beenakker, in Mesoscopic Electron Transport, edited by L. L. 



16 



Sohn, L. P. Kouwenhoven, and G. Schon (Kluwer Academic, Dordrecht, 1997). 

N. S. Bakhvalov, G. S. Kazacha, K. K. Likharev, and S. I. Serdyukova, Physica B 173, 319 

(1991) . 

K. Nabors, S. Kim, and J. White, IEEE Trans, on Microwave Theory and Techniques 40, 1496 

(1992) , software available at ftp://rle-vlsi.mit.edu/pub/fastcap/. 

K. K. Likharev and K. A. Matsuoka, Appl. Phys. Lett. 67, 3037 (1994). 



17 



1.574463e + 07 

1.832335e + 07 e 

1.832340e + 07 . e 

1.832342e + 07 . . e 

1.832412e + 07 . . . e 

1.832451e + 07 . . . . e 

1.832525e + 07 e 

1.832588e + 07 e 

1.832873e + 07 e 

1.833854e + 07 e 

1.834973e + 07 e 

1.844754e + 07 e 

1.861754e + 07 e 

1.899889e + 07 e h 

1.899894e + 07 e h . 

1.899896e + 07 e h . . 

1.899896e + 07 e h . . . 

1.899899e + 07 e h . . . . 

1.899899e + 07 e h 

1.900215e + 07 e h 

1.900377e + 07 e h 

1.902246e + 07 e h 

1.906005e + 07 e h 

1.908176e + 07 e h 

1.918411e + 07 e h 

1.926408e + 07 e h 

1.995202e + 07 e h 

2.001400e + 07 e h 

2.005999e + 07 e h 

2.024492e + 07 e h 

2.025450e + 07 e h 

2.025463e + 07 e h 

2.025538e + 07 e h 

2.025579e + 07 e h 

2.025601e + 07 e h 

2.025622e + 07 e . . . . h 

2.025626e + 07 e . . . h 

2.025627e + 07 e . . h 

2.025629e + 07 e . h 

2.025630e + 07 eh 
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TABLE I. Conduction diagram in a Monte Carlo simulation of an array near a ground plane 

with no background charge, symmetrically biased near threshold, iV = 40, A = 0.96. The first 

column is time, subsequent columns show the charge on each island. 
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3 


4 


5 


6 


7 


8 


9 


10 


11 


12 


1 


2.99525 


0.33124 


0.25056 


0.20843 


0.18008 


0.15821 


0.13997 


0.12453 


0.11055 


0.09774 


2 


0.00456 


0.01429 


0.02383 


0.03347 


0.04319 


0.05314 


0.06333 


0.07441 


0.08573 


0.09768 


3 
4 


2.56553 




















5 


-0.01817 


2.59339 


















6 


-0.00665 


-0.02876 


2.59838 
















7 


-0.00319 


-0.01199 


-0.03172 


2.60029 














8 


-0.00181 


-0.00632 


-0.01391 


-0.03303 


2.60122 












9 


-0.00115 


-0.00385 


-0.00764 


-0.01485 


-0.03372 


2.60174 










10 


-0.00070 


-0.00267 


-0.00473 


-0.00834 


-0.01537 


-0.03412 


2.60205 








11 


-0.00050 


-0.00172 


-0.00323 


-0.00531 


-0.00874 


-0.01567 


-0.03434 


2.60223 






12 


-0.00037 


-0.00122 


-0.00228 


-0.00357 


-0.00562 


-0.00897 


-0.01584 


-0.03448 


2.60231 




13 


-0.00029 


-0.00097 


-0.00163 


-0.00259 


-0.00387 


-0.00576 


-0.00902 


-0.01595 


-0.03452 


2.60231 


14 


-0.00022 


-0.00077 


-0.00130 


-0.00186 


-0.00270 


-0.00405 


-0.00611 


-0.00906 


-0.01591 


-0.03447 


15 


-0.00017 


-0.00055 


-0.00096 


-0.00142 


-0.00200 


-0.00281 


-0.00407 


-0.00585 


-0.00910 


-0.01585 


16 


-0.00013 


-0.00041 


-0.00073 


-0.00107 


-0.00149 


-0.00203 


-0.00276 


-0.00400 


-0.00580 


-0.00897 


17 


-0.00010 


-0.00032 


-0.00055 


-0.00080 


-0.00110 


-0.00149 


-0.00197 


-0.00276 


-0.00382 


-0.00562 


18 


-0.00009 


-0.00026 


-0.00041 


-0.00057 


-0.00076 


-0.00100 


-0.00137 


-0.00190 


-0.00265 


-0.00365 


19 


-0.00002 


-0.00011 


-0.00025 


-0.00042 


-0.00063 


-0.00088 


-0.00121 


-0.00121 


-0.00154 


-0.00223 


20 


-0.00002 


-0.00008 


-0.00015 


-0.00023 


-0.00033 


-0.00046 


-0.00063 


-0.00069 


-0.00092 


-0.00124 


21 


-0.00001 


-0.00003 


-0.00005 


-0.00007 


-0.00010 


-0.00014 


-0.00019 


-0.00021 


-0.00028 


-0.00037 



TABLE II. Partial capacitance matrices generated by FastCap for a cubic array with N = 20 
cubes between large parallel plates. Rows 1 and 2 belong to C, rows 3-21 belong to C. By 
convention, all > 0. All values in 10~ 16 F. The array parameters are b = 1.0, d = 0.04, 
A = 0.356. 
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FIG. 1. ID array of tunnel junctions: (a) general schematic, (b) "ground plane" or "tridiagonal 
model" schematic. 

FIG. 2. Convergence of S/(0)/2e/ calculated by mso as a function of N r , the number of current 
runs in the statistical ensemble, for various values of total transferred charge Q. Results for an 
array near ground plane, N = 20, A = 0.5, at crossover current. 

FIG. 3. Parametric plots of S/(0)/2e/ vs. /. Crossover current I c is defined as the current at 
which Si(0)/2el = 1/s/N. Circles: array near ground plane. Squares: array without ground plane. 
Diamonds and triangles: arrays with random background charge. Dotted lines are only guides for 
the eye. 

FIG. 4. Crossover current I c as a function of the NX, the ratio of array length to charge soliton 
length, in arrays near ground plane with no background charge. 

FIG. 5. T vs. time for a Monte Carlo simulation of an array near ground plane with no back- 
ground charge, symmetrically biased near threshold, N = 20, A = 0.55. Slow points occur before 
the first charge enters the array (a), and when either charge is ~ N/3 junctions away from the 
nearest edge (b,c). Dotted line is only a guide for the eye. 

FIG. 6. Array, without ground plane, between parallel plate externals, (a) Schematics, (b) 
geometric model with surfaces divided into panels for geometric capacitance calculation, and (c) 
closeup of a single island paneling. 

FIG. 7. Crossover current vs. inverse soliton length for arrays without ground plane between 
parallel plate externals. 

FIG. 8. Crossover current as a function of inverse soliton length in arrays near ground plane with 
random background charge. Each point with its error bars represents the average and standard 
deviation over ensembles of fifty different random background charge distributions. Error bar 
widths are scaled for readability: narrowest for N = 10, widest for N = 40. 
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FIG. 9. Typical charge distribution in array with random background charges N = 20, A = 2. 
Filled bars represent %, solid lines represent the sandpile at V = Vt, and dotted lines represent 
conduction charges in a simple alternating e — h scenario. 

FIG. 10. Histograms of the inverse total tunneling rate 1/r at / = I c for an array near ground 
plane (N = 40, A = 0.995) with several different random charge distributions. r q is the seed given 
to the random number generator for background charges. 
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